Il significato sociale della musica in età evolutiva

Laboratorio didattico: esempio di relazione

Motivazione dello studio di Mehr, Song e Spelke (2016)

Mehr, Song e Spelke (2016) discutono i significati sociali della musica in età evolutiva. Secondo questi ricercatori, la musica è capace di trasmettere significati sociali anche nel caso di bambini molto piccoli. Per valutare questa ipotesi, Mehr et al. (2016) esaminano un campione di 32 bambini di età media pari a 5.6 mesi (SD = 0.31, gamma di variazione: 5.1–6.1). Nel primo esperimento, Mehr e al. si chiedono se la familiarità con una melodia possa modulare la preferenza sociale, ovvero i tempi di fissazione dello sguardo (looking-time) che i bambini dirigono verso un individuo adulto sconosciuto che, cantando, produce una melodia familiare oppure non familiare.

Metodo

Nella prima fase dell’esperimento, ai genitori dei bambini veniva chiesto di imparare una nuova ninna-nanna e di ripeterla spesso al figlio/a. Nella fase test dell’esperimento, ai bambini venivano presentati due video, uno di fianco all’altro. Uno dei due video faceva vedere uno sconosciuto che cantava la melodia che, in precedenza, i genitori avevano presentato ai bambini; l’altro video faceva vedere un altro sconosciuto che cantava una ninna-nanna non conosciuta ai bambini. La variabile dipendente era il tempo di fissazione dello sguardo del bambino.

Nel pre-test, ai bambini venivano mostrati gli stessi due adulti utilizzati nella fase test; in questo caso, però, gli attori si limitavano a sorridere.

Risultati

Mediante una prima analisi statistica (relativa ai dati del pre-test) i ricercatori si sono chiesti se i bambini preferissero uno dei due sconosciuti, quando non cantavano. In una seconda analisi statistica, i ricercatori si sono chiesti se i bambini dimostrano una preferenza per lo sconosciuto che canta la ninna-nanna familiare, piuttosto che per lo sconosciuto che canta la ninna-nanna non familiare – ciò dovrebbe corrispondere a tempi di fissazione maggiri per il primo volto rispetto al secondo.

L’analisi statistica eseguita dagli autori è un test \(t\) di Student. L’articolo di Mehr et al. (2016) riporta un valore della statistica test pari a \(t_{31} = 2.96\), un \(p\)-valore di 0.06 e un intervallo di confidenza al 95% pari a [0.529, 0.658]. La dimensione dell’effetto, misurata con il \(d\) di Cohen, è pari a 0.52.

Conclusioni

Mehr e al. (2016) riportano che i bambini di cinque mesi dirigono in modo preferenziale la propria attenzione verso un adulto sconosciuto che canta una melodia familiare piuttosto che verso un adulto sconosciuto che canta una melodia simile, ma sconosciuta. Questi risultati sono stati trovati quando la familiarizzazione con la ninna-nanna target veniva creata mediante l’interazione tra il bambino e un adulto che faceva parte dell’ambiente domestico. L’attenzione dei bambini, invece, non veniva modulata dalla familiarità con la melodia nel caso in cui la melodia in questione fosse stata presentata in precedenza al bambino in un modo diverso, ovvero essendo prodotta da un giocattolo, o da un adulto che il bambino conosceva poco. In base a tali risultati, Mehr e al. (2016) concludono che le melodie prodotte nell’ambiente domestico, nell’interazione tra i bambini e gli adulti, sono dotate di un particolare significato sociale per i bambini.

Scopo del presente progetto

L’analisi dei dati descritta qui di seguito ha lo scopo di replicare alcuni dei risultati riportati da Mehr et al. (2016) utilizzando i dati resi disponibili dagli autori. Le stesse analisi statistiche verranno eseguite usando due metodi diversi: l’approccio frequentista (quello usato nell’articolo) e l’approccio bayesiano.

Analisi dei dati

Per replicare i risultati riportati da Mehr et al. (2016), iniziamo a leggere in R i dati contenuti nel file MehrSongSpelke_exp_1.csv. Questo risultato si ottiene utilizzando la funzione import() contenuta nel pacchetto rio. Utilizziamo anche le funzionalità dei pacchetti tidyverse per manipolare i dati e creare grafici.

Dopo avere cambiato la cartella di lavoro in RStudio, carichiamo i pacchetti che verranno utilizzati mediante le seguenti istruzioni:

suppressPackageStartupMessages({
  library("here")
  library("tidyverse")
  library("patchwork")
  library("bayesplot")
})  

theme_set(bayesplot::theme_default(base_size = 12))
bayesplot::color_scheme_set("brightblue") 

options(
  digits = 3,
  width = 68,
  str = strOptions(strict.width = "cut"),
  crayon.enabled = TRUE
)

knitr::opts_chunk$set(
  comment = "#>",
  collapse = TRUE,
  message = FALSE,
  warning = FALSE,
  error = FALSE,
  width = 68,
  fig.align = "center",
  fig.width = 6,
  fig.asp = 0.618, # 1 / phi
  fig.show = "hold",
  dpi = 300,
  fig.pos = "h",
  cache.extra = knitr::rand_seed,
  tidy.opts = list(width.cutoff = 70),
  tidy = "styler"
)

mehr <- rio::import("MehrSongSpelke_exp_1.csv")

I nomi delle colonne dell’oggetto mehr vengono restituiti usando la funzione names() che prende come argomento il nome del data.frame:

names(mehr)

Consideriamo la variabile exp1:

unique(mehr$exp1)
#> [1] 1 0

Quando exp1 assume il valore 1, i dati appartengono al primo esperimento, quello di interesse qui. Per selezionare i dati del primo esperimento usiamo la funzione filter():

df_exp1 <- mehr %>% 
  dplyr::filter(exp1 == 1)

L’istruzione precedente dice che l’oggetto mehr (che è un DataFrame) viene passato alla funzione filter() contenuta nel pacchetto dplyr. L’argomento passato a filter() istruisce R a selezionare solo le righe del data.frame mehr nelle quali la colonna exp ha valore 1. Il sottoinsieme delle righe selezionate corrispone alle osservazioni che fanno parte del primo esperimento. L’operatore di attribuzione <- assegna al risultato che abbiamo ottenuto (ovvero, ai dati del primo esperimento) il nome df_exp1.

La variabile Test_Proportion_Gaze_to_Singer corrisponde al tempo di fissazione dello sguardo (espresso nei termini di una proporzione). Costruiamo un istogramma:

tibble(y = df_exp1$Test_Proportion_Gaze_to_Singer) %>% 
ggplot(aes(y)) +
  geom_histogram(
    aes(y = ..density..)
  )

Più utile è una stima della densità della frequenza dei dati:

tibble(y = df_exp1$Test_Proportion_Gaze_to_Singer) %>% 
ggplot(aes(y)) +
  geom_density()

Media e deviazione standard delle proporzioni di tempo trascorso a fissare il volto target (sconosciuto che canta la ninna-nanna familiare) si ottengono nel modo seguente:

df_exp1 %>% 
  summarise(
    p = mean(Test_Proportion_Gaze_to_Singer),
    std = sd(Test_Proportion_Gaze_to_Singer)
  )
#>       p   std
#> 1 0.593 0.179

Approccio frequentista

Test \(t\) di Student per un solo campione

Utilizziamo ora l’approccio frequentista per fare inferenza su \(\mu\), ovvero sulla proporzione media del tempo di fissazione rivolto allo sconosciuto target, piuttosto che verso lo sconosciuto che canta una ninna-nanna ignota al bambino.

L’ipotesi nulla è

\[ H_0: \mu = 0.5. \]

Se non vi è alcuna preferenza, la variabile Test_Proportion_Gaze_to_Singer dovrebbe assumere il valore 0.5. Se invece i bambini preferiscono la melodia familiare, il valore medio di Test_Proportion_Gaze_to_Singer dovrebbe essere maggiore di 0.5; se i bambini preferiscono la melodia sconosciuta, il valore medio di Test_Proportion_Gaze_to_Singer dovrebbe essere minore di 0.5.

Iniziamo ad esaminare la distribuzione dei dati. Usiamo il diagramma quantile-quantile per verificare l’ipotesi di gaussianità (necessaria per eseguire il test statistico riportato dagli autori):

tibble(y = df_exp1$Test_Proportion_Gaze_to_Singer) %>% 
  ggplot(aes(sample = y)) +
  stat_qq() + 
  stat_qq_line()

Dato che l’assunzione di gaussianità è ragionevole, procediamo eseguendo il test \(t\) di Student. La proporzione del looking-time diretta allo sconosciuto che canta la melodia familiare è

m <- mean(df_exp1$Test_Proportion_Gaze_to_Singer)
m
#> [1] 0.593

pari a 0.59. Abiamo trovato tale valore usando la funzione mean() che prende, come argomento, il vettore che contiene la proporzione del looking-time diretta allo sconosciuto che canta la melodia familiare per ciascun bambino nel primo esperimento. Utilizziamo la sintassi $ per estrarre dal DataFrame df_exp1 la colonna di dati Test_Proportion_Gaze_to_Singer.

L’errore standard della media è

\[ \sigma_{\bar{y}} = \frac{\sigma}{\sqrt{n}} \] ovvero

se <- sd(df_exp1$Test_Proportion_Gaze_to_Singer) / sqrt(length(df_exp1$Test_Proportion_Gaze_to_Singer))
se
#> [1] 0.0316

La funzione length() ritorna il numero di elementi di un vettore ovvero, nel nostro caso, l’ampiezza del campione.

Applichiamo ora la formula per il calcolo della statistica T, ovvero

\[ T = \frac{\bar{y} - \mu_0}{\hat{\sigma}_{\bar{y}}} = \frac{0.59 - 0.50}{0.032} = 2.96 \]

ovvero

T <- (m - 0.5) / se
T
#> [1] 2.96

I gradi di libertà sono:

dof <- length(df_exp1$Test_Proportion_Gaze_to_Singer) - 1
dof
#> [1] 31

Per un test bidirezionale, il \(p\)-valore diventa

round(2 * (1 - pt(T, dof)), 3)
#> [1] 0.006

Nell’istruzione precedente abbiamo calcolato l’area sottesa alla ditribuzione \(t_{31}\) nell’intervallo \([2.960, \infty]\) e l’abbiamo moltiplicata per 2. Per trovare l’area nell’intervallo \([2.960, \infty]\), abbiamo calcolato il valore della funzione di ripartizione in corrispondenza del valore 2.960 (ovvero l’area nell’intervallo \([-\infty, 2.960]\)) e abbiamo sottratto tale valore da 1.0.

Gli stessi risultati, che coincidono con quelli trovati sopra e con quelli riportati dagli autori, si possono facilmente ottenere usando la funzione t.test()

t.test(df_exp1$Test_Proportion_Gaze_to_Singer, mu = 0.5)
#> 
#>  One Sample t-test
#> 
#> data:  df_exp1$Test_Proportion_Gaze_to_Singer
#> t = 3, df = 31, p-value = 0.006
#> alternative hypothesis: true mean is not equal to 0.5
#> 95 percent confidence interval:
#>  0.529 0.658
#> sample estimates:
#> mean of x 
#>     0.593

L’intervallo di confidenza al 95% si calcola come

\[ \text{stima del parametro} \pm t \times \text{errore standard} \] Nel caso presente, abbiamo:

round(m + c(-1, 1) * qt(0.975, dof) * se, 3)
#> [1] 0.529 0.658

Tale risultato riproduce quello riportato da Mehr et al. (2016).

Una stima della dimensione dell’effetto si ottiene come il rapporto tra la differenza tra le medie (in questo caso, proporzioni) e la relativa deviazione standard:

round(
  (m - 0.5) / sd(df_exp1$Test_Proportion_Gaze_to_Singer),
  3
)
#> [1] 0.523

Anche questo risultato coincide con il valore riportato dagli autori.

Analisi Bayesiana

Ci poniamo ora il problema di ripetere le analisi statiche precedenti utilizzando l’approccio bayesiano. Carichiamo i pacchetti necessari:

suppressPackageStartupMessages({
  library("brms")
  library("rstan")
  library("cmdstanr")
  library("posterior")
  library("loo")
})

rstan_options(auto_write = TRUE) # avoid recompilation of models
options(mc.cores = parallel::detectCores()) # parallelize across all CPUs
Sys.setenv(LOCAL_CPPFLAGS = "-march=native") # improve execution time
rstan_options(auto_write = TRUE) # avoid recompilation of models
options(mc.cores = parallel::detectCores()) # parallelize across all CPUs
Sys.setenv(LOCAL_CPPFLAGS = "-march=native") # improve execution time

I dati sono i seguenti:

y <- sort(df_exp1$Test_Proportion_Gaze_to_Singer)
y
#>  [1] 0.263 0.282 0.366 0.418 0.419 0.437 0.460 0.462 0.473 0.499
#> [11] 0.500 0.508 0.508 0.531 0.542 0.542 0.572 0.586 0.601 0.603
#> [21] 0.624 0.683 0.700 0.724 0.760 0.763 0.777 0.789 0.811 0.900
#> [31] 0.938 0.951

Considereremo \(y\) (così come hanno fatto gli autori) come una variabile continua per la quale è sensato assumere un meccanismo generatore dei dati gaussiano di media e deviazione standard ignote. L’inferenza riguarda il parametro \(\mu\) di tale distribuzione gaussiana.

Funzione brm()

Utilizzando delle distribuzioni a priori per i parametri debolmente informative, facciamo inferenza utilizzando la funzione brm() del pacchetto brms:

fit <- brm(
  data = df_exp1, 
  family = gaussian(),
  Test_Proportion_Gaze_to_Singer ~ 1,
  prior = c(
    prior(normal(0, 2.0), class = Intercept),
    prior(cauchy(0, 2.0), class = sigma)
  ),
  iter = 4000, 
  refresh = 0, 
  chains = 4,
  backend = "cmdstanr"
)
#> Running MCMC with 4 chains, at most 8 in parallel...
#> 
#> Chain 1 finished in 0.1 seconds.
#> Chain 2 finished in 0.1 seconds.
#> Chain 3 finished in 0.1 seconds.
#> Chain 4 finished in 0.1 seconds.
#> 
#> All 4 chains finished successfully.
#> Mean chain execution time: 0.1 seconds.
#> Total execution time: 0.3 seconds.

Le stime a posteriori dei parametri si recuperano nel modo seguente:

fixef(fit)
#>           Estimate Est.Error  Q2.5 Q97.5
#> Intercept    0.594    0.0336 0.527 0.659

Nell’output, Estimate corrisponde alla stima della media a posteriori, mentre Q2.5 e Q97.5 forniscono gli estremi dell’intervallo di credibilità al 95%. I risultati sono simili a quelli riportati dagli autori.

Modello in linguaggio Stan

Ripetiamo ora l’analisi precedente utilizzano un modello scritto in linguaggio Stan. Il modello riportato di seguito verrà salvato nella cartella code con il nome normalmodel.stan:

modelString = "
data {
  int<lower=0> N;
  vector[N] y;
}
parameters {
  real mu;
  real<lower=0> sigma;
}
model {
  mu ~ normal(0.5, 0.5);
  sigma ~ cauchy(0, 1);
  y ~ normal(mu, sigma);
}
"
writeLines(modelString, con = "code/normalmodel.stan")

Sistemiamo i dati nel formato appropriato per Stan:

data_list <- list(
  N = length(df_exp1$Test_Proportion_Gaze_to_Singer),
  y = df_exp1$Test_Proportion_Gaze_to_Singer
)
data_list
#> $N
#> [1] 32
#> 
#> $y
#>  [1] 0.603 0.683 0.724 0.282 0.499 0.951 0.418 0.938 0.500 0.586
#> [11] 0.473 0.508 0.811 0.572 0.777 0.263 0.508 0.437 0.542 0.601
#> [21] 0.419 0.789 0.760 0.624 0.366 0.462 0.900 0.531 0.542 0.700
#> [31] 0.763 0.460

Leggiamo il file in cui abbiamo salvato il codice Stan

file <- file.path("code", "normalmodel.stan")

compiliamo il modello

mod <- cmdstan_model(file)

ed eseguiamo il campionamento MCMC:

fit <- mod$sample(
  data = data_list,
  iter_sampling = 4000L,
  iter_warmup = 2000L,
  seed = 123,
  chains = 4L,
  parallel_chains = 2L,
  refresh = 0,
  thin = 1
)

La convergenza e il “mixing” del campionamanto MCMC possono essere controllate mediante il trace plot che mostra l’andamento delle simulazioni e ci dice se stiamo effettivamente utilizzando una distribuzione limite:

Convertiamo l’oggetto fit in formato stanfit,

stanfit <- rstan::read_stan_csv(fit$output_files())

Esaminiamo i trace-plot:

stanfit %>% 
  mcmc_trace(pars = c("mu", "sigma"), size = 0.1)

La stima della distribuzione a posteriori per \(\mu\) per ciascuna delle quattro catene usate può essere rappresentata graficamente nel modo seguente:

mcmc_dens_overlay(stanfit, pars = "mu") + 
  ylab("density")

Esaminiamo la statistica \(\hat{R}\):

bayesplot::rhat(stanfit, pars = c("mu", "sigma"))
#>    mu sigma 
#>     1     1

Utilizzando l’oggetto stanfit, possiamo recuperare la statistica di Geweke nel modo seguente:

fit_mcmc <- As.mcmc.list(
  stanfit,
  pars = c("mu", "sigma")
)
coda::geweke.diag(fit_mcmc, frac1 = .1, frac2 = .5) 
#> [[1]]
#> 
#> Fraction in 1st window = 0.1
#> Fraction in 2nd window = 0.5 
#> 
#>    mu sigma 
#>  1.21  1.60 
#> 
#> 
#> [[2]]
#> 
#> Fraction in 1st window = 0.1
#> Fraction in 2nd window = 0.5 
#> 
#>    mu sigma 
#> 0.631 0.874 
#> 
#> 
#> [[3]]
#> 
#> Fraction in 1st window = 0.1
#> Fraction in 2nd window = 0.5 
#> 
#>     mu  sigma 
#> -0.271 -0.589 
#> 
#> 
#> [[4]]
#> 
#> Fraction in 1st window = 0.1
#> Fraction in 2nd window = 0.5 
#> 
#>     mu  sigma 
#> -0.276  0.180

La statistica di Geweke è uguale a zero quando le medie delle due porzioni della catena di Markov sono uguali. Valori maggiori di \(\mid 2 \mid\) suggeriscono che la catena non ha ancora raggiunto una distribuzione stazionaria.

Le stime a posteriori dei parametri si ottengono con

s <- summary(
  stanfit, 
  c("mu", "sigma"), 
  probs = c(0.025, 0.50, 0.975),
  use_cache = TRUE
  )
s$summary  # all chaines merged
#>        mean  se_mean     sd  2.5%   50% 97.5% n_eff Rhat
#> mu    0.594 0.000309 0.0332 0.529 0.594 0.659 11504    1
#> sigma 0.186 0.000235 0.0251 0.144 0.183 0.241 11427    1

I risultati replicano quelli ottenuti in precedenza.

Modifichiamo ora il modello inserendo una distribuzione a priori più sensata per il parametro \(\mu\). Essendo una proporzione, deve essere un numero compreso tra 0 e 1. Pertanto, quale distribuzione a priori, è sensato usare una distribuzione Beta. Utilizzeremo qui una Beta(2,2), ovvero una distribuzione a priori debolmente informativa.

modelString = "
data {
  int<lower=0> N;
  vector[N] y;
}
parameters {
  real mu;
  real<lower=0> sigma;
}
model {
  mu ~ beta(2, 2);
  sigma ~ cauchy(0, 1);
  y ~ normal(mu, sigma);
}
"
writeLines(modelString, con = "code/normalmodel2.stan")

Leggiamo il file in cui abbiamo salvato il codice Stan

file2 <- file.path("code", "normalmodel2.stan")

compiliamo il modello

mod2 <- cmdstan_model(file2)

ed eseguiamo il campionamento MCMC:

fit2 <- mod2$sample(
  data = data_list,
  iter_sampling = 4000L,
  iter_warmup = 2000L,
  seed = 123,
  chains = 4L,
  parallel_chains = 2L,
  refresh = 0,
  thin = 1
)

Troviamo le stime a posteriori dei parametri:

stanfit2 <- rstan::read_stan_csv(fit2$output_files())

s <- summary(
  stanfit2, 
  c("mu", "sigma"), 
  probs = c(0.025, 0.50, 0.975),
  use_cache = TRUE
  )
s$summary  # all chaines merged
#>        mean  se_mean     sd  2.5%   50% 97.5% n_eff Rhat
#> mu    0.593 0.000293 0.0334 0.527 0.594 0.659 13016    1
#> sigma 0.186 0.000227 0.0250 0.145 0.183 0.242 12073    1

I risultati sono quasi identici a quelli trovati in precedenza.

Possiamo dunque concludere, con un grado di certezza soggettiva del 95%, che siamo sicuri che la proporzione del tempo complessivo di fissazione dello sguardo verso lo sconosciuto che cantava la melodia familiare (piuttosto che verso lo sconosciuto che cantava una ninn-nanna ignota), è compresa nell’intervallo [0.527, 0.659].

Confronto tra due gruppi

Approccio frequentista (test \(t\) di Student per campioni appaiati)

In una diversa analisi statistica, Mehr e al. (2016) hanno confrontato il tempo di fissazione dello sguardo diretto verso lo sconosciuto che sorride soltanto e verso lo sconosciuto che canta la ninna-nanna familiare. Ci dovremmo aspettare un tempo di fissazione maggiore nel secondo caso rispetto al primo.

La distribuzione dei tempi di fissazione dello sguardo nelle due condizioni è presentata nel pannello a della Figura 2 del loro articolo. Riproduciamo qui sotto la figura. I tempi di fissazione dello sguardo nelle condizioni baseline e test corrispondono, rispettivamente, alle variabili Baseline_Proportion_Gaze_to_Singer e Test_Proportion_Gaze_to_Singer.

Per replicare la Figura 2a riportata nell’articolo dobbiamo sistemare i dati nel formato ``long’’, ovvero nel formato in cui a ciascuna colonna del DataFrame corrisponde una variabile. Nel caso presente, possiamo creare una variabile chiamata condizione, con modalità baseline e test, e una seconda variabile, y, che corrisponde ai tempi di fissazione dello sguardo. Per creare un DataFrame che contiene queste due variabili, procediamo come segue.

condizione <- factor(
  rep(
    c("baseline", "test"), 
    each = length(df_exp1$Baseline_Proportion_Gaze_to_Singer)
    )
  )
y <- c(df_exp1$Baseline_Proportion_Gaze_to_Singer, df_exp1$Test_Proportion_Gaze_to_Singer)
df2 <- data.frame(condizione, y)
head(df2)
#>   condizione     y
#> 1   baseline 0.437
#> 2   baseline 0.413
#> 3   baseline 0.754
#> 4   baseline 0.439
#> 5   baseline 0.475
#> 6   baseline 0.871
dim(df2)
#> [1] 64  2

Il diagramma a scatola (simile a quello riportato nell’articolo) si crea con le seguenti istruzioni:

p <- df2 %>%
  ggplot(aes(x = condizione, y = y, fill = condizione)) +
  geom_boxplot(alpha = 0.7) +
  scale_y_continuous(
    name = "Proporzione di fissazioni verso il \ncantante di canzone familiare"
  ) +
  scale_x_discrete(name = "Condizione") +
  theme(legend.position = "none") +
  geom_hline(yintercept = 0.5, linetype = "dashed") +
  stat_summary(
    fun.y = mean, colour = "darkred", geom = "point",
    shape = 18, size = 6
  )
p

Una rappresentazione più efficace dei dati si ottiene riportando, per ciascuna condizione, tutte le singole osservazioni, insieme alla mediana di ciascun gruppo:

 mehr_summary <- df2 %>%
  group_by(condizione) %>%
  summarize(
    lt_mean = mean(y),
    lt_sd = sd(y),
    lt_median = median(y)
  ) %>%
  ungroup()

df2 %>%
  ggplot(
    aes(x = condizione, y = y, color = condizione)
  ) +
  ggforce::geom_sina(aes(color = condizione, size = 3, alpha = .5)) +
  geom_errorbar(
    aes(y = lt_median, ymin = lt_median, ymax = lt_median),
    data = mehr_summary, width = 0.5, size = 3
  ) +
  labs(
    x = "",
    y = "Proporzione di fissazioni verso il \ncantante di canzone familiare"
  ) +
  theme(legend.position = "none") +
  scale_colour_grey(start = 0.7, end = 0)

Un confronto tra queste due condizioni viene effettuato con un test \(t\) di Student per campioni appaiati. Ovvero

diff <- df_exp1$Test_Proportion_Gaze_to_Singer - df_exp1$Baseline_Proportion_Gaze_to_Singer
diff
#>  [1]  0.16561  0.27049 -0.03035 -0.15722  0.02390  0.08002  0.18104
#>  [8]  0.17894  0.08367 -0.21324  0.09395 -0.18951  0.21778 -0.04311
#> [15]  0.16254 -0.05399  0.19752 -0.06739  0.07277  0.09682 -0.14536
#> [22]  0.53284  0.06011  0.24154 -0.00545  0.17707  0.13171  0.05731
#> [29] -0.27932  0.11023  0.34093  0.02479

Verifichiamo la gaussianità dei dati:

tibble(y = diff) %>% 
  ggplot(aes(sample = y)) +
  stat_qq() + 
  stat_qq_line()

La media delle differenze è:

m <- mean(diff)
m
#> [1] 0.0724

L’errore standard delle differenze è:

se <- sd(diff) / sqrt(length(diff))
se
#> [1] 0.03

La statistica \(T\) è:

T <- (m - 0) / se
T
#> [1] 2.42

I gradi di libertà sono:

nu <- length(diff) - 1
nu
#> [1] 31

Il p-valore è:

2 * (1 - pt(T, nu))
#> [1] 0.0218

Usando t.test() otteniamo:

t.test(diff, mu = 0)
#> 
#>  One Sample t-test
#> 
#> data:  diff
#> t = 2, df = 31, p-value = 0.02
#> alternative hypothesis: true mean is not equal to 0
#> 95 percent confidence interval:
#>  0.0113 0.1335
#> sample estimates:
#> mean of x 
#>    0.0724

Questi risultati replicano i valori riportati dagli autori: \(t_{31}\) = 2.42, \(p\) = 0.22.

L’intervallo di confidenza è

mean(diff) + c(-1, 1) * qt(.975, nu) * se
#> [1] 0.0113 0.1335

La dimensione dell’effetto è

d <- mean(diff) / sd(diff)
d
#> [1] 0.427

Le stime dell’intervallo di confidenza e della dimensione dell’effetto coincidono con i valori riportati da Mehr e al. (2016).

Analisi bayesiana

Dato che, nel caso di campioni appaiati, abbiamo un unico valore della varibile dipendente per ciascuna “coppia”, l’analisi statistica è sostanzialmente identica a quella già descritta sopra. In questo caso, la variabile dipendente è chiamata diff.

data2_list <- list(
  N = length(diff),
  y = diff
)
data2_list
#> $N
#> [1] 32
#> 
#> $y
#>  [1]  0.16561  0.27049 -0.03035 -0.15722  0.02390  0.08002  0.18104
#>  [8]  0.17894  0.08367 -0.21324  0.09395 -0.18951  0.21778 -0.04311
#> [15]  0.16254 -0.05399  0.19752 -0.06739  0.07277  0.09682 -0.14536
#> [22]  0.53284  0.06011  0.24154 -0.00545  0.17707  0.13171  0.05731
#> [29] -0.27932  0.11023  0.34093  0.02479

Usiamo nuovamente il modello normalmodel2.stan che abbiamo già compilato in precedenza ed eseguiamo il campionamento MCMC:

fit3 <- mod2$sample(
  data = data2_list,
  iter_sampling = 4000L,
  iter_warmup = 2000L,
  seed = 123,
  chains = 4L,
  parallel_chains = 2L,
  refresh = 0,
  thin = 1
)

Le stime a posteriori dei parametri si ottengono con:

stanfit3 <- rstan::read_stan_csv(fit3$output_files())

s <- summary(
  stanfit3, 
  c("mu", "sigma"), 
  probs = c(0.025, 0.50, 0.975),
  use_cache = TRUE
  )
s$summary  # all chaines merged
#>        mean  se_mean     sd   2.5%   50% 97.5% n_eff Rhat
#> mu    0.085 0.000281 0.0286 0.0322 0.084 0.146 10406    1
#> sigma 0.176 0.000233 0.0238 0.1374 0.174 0.231 10392    1

I risultati sono molto simili ai precedenti.

Dall’output precedente possiamo ottenere una stima della dimensione dell’effetto:

s$summary[1, 1] /   s$summary[2, 1]
#> [1] 0.482

Anche in questo caso, il risultato è simile a quello ottenuto sopra.

Conclusioni

Utilizzando i dati grezzi messi a disposizione dagli autori siamo stati in grado di replicare i risultati riportati dall’articolo di Mehr e al. (2016). In questa relazione ci siamo focalizzati su due problemi: l’inferenza su una media e l’inferenza su due medie nel caso di dati appaiati. Nel caso presente, i risultati dell’analisi bayesiana sono praticamente identici a quelli prodotti dall’approccio frequentista, se vengono usate distribuzioni a priori debolmente informative. Utilizzando ggplot() siamo anche stati in grado di replicare una delle figure riportate nell’articolo di Mehr e al. (2016). Abbiamo anche proposto una nuova versione di tale figura che sembra più informativa di quella proposta dagli autori.

Le analisi statistiche descritte da Mehr e al. (2016) sono convincenti: i dati non mostrano anomalie degne di nota e le assunzioni alla base dei test statistici svolti dagli autori – ovvero, la gaussianità dei dati – sono sensate. È dunque sensato concludere che la musica sembra possedere un importante valore sociale anche nel caso di bambini molto piccoli.

Riferimenti bibliografici

Mehr, S. A., Song, L. A., & Spelke, E. S. (2016). For 5-month-old infants, melodies are social. Psychological Science, 27(4), 486-501.

Informazioni sulla sessione di lavoro

sessionInfo()
#> R version 4.1.2 (2021-11-01)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Big Sur 10.16
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
#> 
#> locale:
#> [1] it_IT.UTF-8/it_IT.UTF-8/it_IT.UTF-8/C/it_IT.UTF-8/it_IT.UTF-8
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods  
#> [7] base     
#> 
#> other attached packages:
#>  [1] loo_2.4.1            posterior_1.2.0      cmdstanr_0.4.0.9000 
#>  [4] rstan_2.21.3         StanHeaders_2.21.0-7 brms_2.16.3         
#>  [7] Rcpp_1.0.8           bayesplot_1.8.1      patchwork_1.1.1     
#> [10] forcats_0.5.1        stringr_1.4.0        dplyr_1.0.7         
#> [13] purrr_0.3.4          readr_2.1.1          tidyr_1.1.4         
#> [16] tibble_3.1.6         ggplot2_3.3.5        tidyverse_1.3.1     
#> [19] here_1.0.1          
#> 
#> loaded via a namespace (and not attached):
#>   [1] readxl_1.3.1         backports_1.4.1      plyr_1.8.6          
#>   [4] igraph_1.2.11        splines_4.1.2        crosstalk_1.2.0     
#>   [7] TH.data_1.1-0        rstantools_2.1.1     inline_0.3.19       
#>  [10] digest_0.6.29        htmltools_0.5.2      rsconnect_0.8.25    
#>  [13] fansi_1.0.2          magrittr_2.0.2       checkmate_2.0.0     
#>  [16] tzdb_0.2.0           openxlsx_4.2.5       modelr_0.1.8        
#>  [19] RcppParallel_5.1.5   matrixStats_0.61.0   R.utils_2.11.0      
#>  [22] sandwich_3.0-1       xts_0.12.1           prettyunits_1.1.1   
#>  [25] colorspace_2.0-2     rvest_1.0.2          haven_2.4.3         
#>  [28] xfun_0.29            callr_3.7.0          crayon_1.4.2        
#>  [31] jsonlite_1.7.3       lme4_1.1-27.1        survival_3.2-13     
#>  [34] zoo_1.8-9            glue_1.6.1           polyclip_1.10-0     
#>  [37] gtable_0.3.0         emmeans_1.7.2        distributional_0.3.0
#>  [40] R.cache_0.15.0       pkgbuild_1.3.1       abind_1.4-5         
#>  [43] scales_1.1.1         mvtnorm_1.1-3        DBI_1.1.2           
#>  [46] miniUI_0.1.1.1       xtable_1.8-4         foreign_0.8-82      
#>  [49] stats4_4.1.2         DT_0.20              htmlwidgets_1.5.4   
#>  [52] httr_1.4.2           threejs_0.3.3        ellipsis_0.3.2      
#>  [55] pkgconfig_2.0.3      R.methodsS3_1.8.1    farver_2.1.0        
#>  [58] sass_0.4.0           dbplyr_2.1.1         utf8_1.2.2          
#>  [61] tidyselect_1.1.1     labeling_0.4.2       rlang_1.0.0.9000    
#>  [64] reshape2_1.4.4       later_1.3.0          munsell_0.5.0       
#>  [67] cellranger_1.1.0     tools_4.1.2          cli_3.1.1           
#>  [70] generics_0.1.1       broom_0.7.12         ggridges_0.5.3      
#>  [73] evaluate_0.14        fastmap_1.1.0        yaml_2.2.2          
#>  [76] processx_3.5.2       knitr_1.37           fs_1.5.2            
#>  [79] zip_2.2.0            nlme_3.1-155         projpred_2.0.2      
#>  [82] mime_0.12            R.oo_1.24.0          xml2_1.3.3          
#>  [85] compiler_4.1.2       shinythemes_1.2.0    rstudioapi_0.13     
#>  [88] gamm4_0.2-6          curl_4.3.2           reprex_2.0.1        
#>  [91] tweenr_1.0.2         bslib_0.3.1          stringi_1.7.6       
#>  [94] highr_0.9            ps_1.6.0             Brobdingnag_1.2-6   
#>  [97] lattice_0.20-45      Matrix_1.4-0         nloptr_2.0.0        
#> [100] styler_1.6.2         markdown_1.1         shinyjs_2.1.0       
#> [103] tensorA_0.36.2       vctrs_0.3.8          pillar_1.6.5        
#> [106] lifecycle_1.0.1      jquerylib_0.1.4      bridgesampling_1.1-2
#> [109] estimability_1.3     data.table_1.14.2    httpuv_1.6.5        
#> [112] R6_2.5.1             promises_1.2.0.1     gridExtra_2.3       
#> [115] rio_0.5.29           codetools_0.2-18     boot_1.3-28         
#> [118] MASS_7.3-55          colourpicker_1.1.1   gtools_3.9.2        
#> [121] assertthat_0.2.1     rprojroot_2.0.2      withr_2.4.3         
#> [124] shinystan_2.5.0      multcomp_1.4-18      mgcv_1.8-38         
#> [127] parallel_4.1.2       hms_1.1.1            grid_4.1.2          
#> [130] minqa_1.2.4          coda_0.19-4          rmarkdown_2.11      
#> [133] ggforce_0.3.3        shiny_1.7.1          lubridate_1.8.0     
#> [136] base64enc_0.1-3      dygraphs_1.1.1.6